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Abstract 

In  this  paper  we  propose  a  variational  approach  to  a  path  planning 
problem  in  2  dimensions  using  a  level  set  framework.  After  defining  an 
energy  integral  over  the  path,  we  use  gradient  flow  on  the  defined  energy 
and  evolve  the  entire  path  until  a  locally  optimal  steady  state  is  reached. 
Unlike  typical  level  set  implementations  where  the  interface  being  tracked 
is  a  codimension- 1  set,  we  allow  for  paths  with  positive,  varying  widths. 
Applications  of  this  method  extend  to  robotic  motion,  tool-path  milling, 
and  arial  search  patterns  for  example.  Numerical  methods  and  algorithms 
are  given,  and  examples  are  presented. 


1  Introduction 

The  general  problem  of  finding  the  optimal  path  through  a  domain  under  some 
given  constraints  has  multiple  applications.  From  path  planning  for  autonomous 
vehicles  [23],  to  computing  tool  path  trajectories  [8],  to  maximizing  visibility  [5], 
there  exist  wide  applications  for  the  general  solutions  to  this  problem.  Given 
that  the  domain  is  not  homogeneous,  i.e.  there  is  an  associated  cost  function 
to  the  path  in  the  domain,  the  general  solution  begins  to  increase  rapidly  in 
complexity. 

A  specific  instance  of  the  general  problem  is  that  of  finding  an  optimal- 
path  map  for  a  known  environment.  The  optimal-path  map  for  a  known  two- 
dimensional  terrain  is  a  function  u(x,y)  whose  values  describe  how  to  best 
reach  the  goal  from  the  location  ( x ,  y).  Optimality  in  this  case  could  be  shortest 
path,  least  visibility  from  above,  largest  patrol  area,  etc.  Previous  work  on  true 
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optimal-path  maps  for  autonomous  robotics  have  almost  always  investigated 
restricted  cases.  While  a  thorough  search  is  performed,  it  is  usually  over  a  finite 
set  of  locations.  Searching  over  an  infinite  set  of  locations  does  not  allow  the 
possibility  of  reducing  the  problem  to  graph  search. 

Robotic  Motion  In  [23],  the  problem  of  constructing  optimal-path  maps 
for  two-dimensional  polygonal  weighted-region  domains  was  investigated  for 
piecewise-constant  solutions.  The  complexity  of  the  algorithm  rendered  it  in¬ 
compatible  for  real-time  path  generation,  but  when  a  large  amount  of  planning 
time  is  available  (daily  routes  of  a  sentry,  etc.)  it  became  desirable.  In  [15]  path 
planning  for  robots  was  studied  using  level  sets  where  there  were  objects  to  be 
avoided  in  the  domain.  The  method  of  solution  was  to  construct  a  weighted 
distance  function  over  the  entire  domain  and  then,  from  a  final  position,  back 
propogate  the  solution  perpendicular  to  the  level  sets  of  the  distance  function, 
resulting  in  an  optimally  shortest  path.  Path  planning  algorithms  for  mobile 
robots  are  also  described  in  [17], [2], [16].  Also,  in  the  context  of  manipulators 
there  has  been  path  planning  research  done  within  a  variational  framework  [24] . 

Tool  Paths  Other  instances  of  the  general  problem  arise  in  the  generation  of 
tool  paths  for  pocket  machining.  Here,  material  is  milled  from  a  blank,  layer  by 
layer,  until  the  object  is  machined  into  a  manufactured  part.  The  tool  path  for 
a  layer  of  a  pocket  is  the  centerline  path  along  which  a  tool  cuts  the  material. 
Optimal  for  this  problem  involves  minimal  length  and  minimal  curvature  of 
the  path  (high  curvature  induces  more  wear  and  tear  on  the  drill).  In  [1],  an 
approximate  solution  to  the  problem  of  how  to  best  produce  tool  paths  for 
pocket  machining  was  presented.  There,  the  authors  solved  an  elliptic  partial 
differential  equation  (PDE)  boundary  value  problem,  and  used  the  contours  of 
that  solution  to  construct  a  solution.  In  [8]  an  extensive,  categorized,  reference 
list  is  presented  of  papers  related  to  numerical  control  tool  path  generation. 
Some  of  the  areas  of  research  are  isoparametric  paths,  non-isoparametric  paths, 
planar  pocketing  paths,  roughing  paths,  and  space- filling  curve  based  tool  paths. 

Dynamic  Visibility  In  [27]  the  framework  for  studying  visibility  and  its  dy¬ 
namics  using  level  sets  was  established.  This  included  the  use  of  implicit  surfaces 
to  represent  visible  regions  and  the  introduction  of  PDEs  that  govern  horizon 
and  boundary  curves.  In  [5]  various  variational  problems  were  approached  using 
the  framework  established  in  [27],  including  single  and  multiple  point  visibility 
optimization,  the  effects  of  weighted  spatial  regions,  and  the  effects  of  memory. 
In  [5]  a  parameterized  path  planning  algorithm  was  introduced  that  treats  the 
path  as  a  finite  union  of  multiple  observers  which  are  evolved  so  as  to  maximize 
the  accumulated  visibility  along  the  path.  See  also  [29]  for  a  path  planning 
algorithm  based  on  visibility. 

Related  ideas  in  image  processing  having  to  do  with  active  contours  and 
snakes  can  be  found  [14], [3], [11].  These  papers  introduce  curve  evolutions  re¬ 
sulting  from  variational  frameworks  used  to  segment  images.  Often  these  imple- 
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mentations  involve  the  maximization  of  a  quantity  defined  by  the  image,  such 
as  the  norm  of  the  gradient,  over  a  curve,  and  they  often  involve  regularizing 
terms  resulting  in  mean  curvature  flow  when  in  a  level  set  setting. 

Something  the  above  examples  have  in  common  that  it  is  not  only  the  path 
that  is  optimized  but  also  the  set  that  the  path  carves  out  that  must  be  con¬ 
sidered.  For  the  robot  path  planning  algorithm  the  width  of  the  set  may  be 
dependent  upon  the  size  of  the  robot.  In  the  tool  planning  problem  it  can  be 
the  width  of  the  drill  bit.  For  searching  we  can  consider  either  the  region  that 
the  observer  can  see  to  be  analogous  to  the  size  of  the  drill  bit  in  the  milling 
example.  A  reverse  case  would  be  if  an  eluder  did  not  want  to  be  seen,  and  then 
the  size  of  the  region  where  the  eluder  can  be  detected  would  be  analogous  to 
the  drill  bit  size. 

In  this  paper,  we  investigate  the  general  problem  of  finding  a  “search  path” 
through  a  domain  where  we  know  some  information  about  where  targets  may 
be  located.  An  agent  searching  such  a  domain  would  want  to  have  a  path  that 
satisfies  being  shortest  with  having  a  high  confidence  of  finding  targets.  The 
searcher  can  only  “see”  a  finite  distance  about  it  at  any  given  point,  and  this 
distance  may  vary  spatially  according  to  local  weather  conditions,  altitude,  etc. 
Therefore,  we  wish  to  generate  an  optimal  path  that  gives  a  certain  level  of 
confidence  of  locating  targets  which  will  be  determined  via  the  information  we 
know  about  the  domain. 

The  remainder  of  the  paper  is  organized  as  follows,  in  the  next  section, 
we  formulate  the  search  path  problem  in  a  general  framework,  with  general 
metrics  describing  the  optimization,  and  constraints  stemming  from  the  search 
path  problem.  Following  that,  we  introduce  the  level  set  method,  and  then  our 
algorithm.  We  present  simulations  of  canonical  examples  demonstrating  the 
method  and  conclude  with  some  remarks  about  the  generality  of  the  method. 


2  Problem  Formulation 


The  general  path  planning  problem  has  had  many  formulations.  For  a  given  set 
17  £  R",  we  seek  a  path  rsE,  with  the  following  properties: 

1.  Optimize  some  function  of  T  (Arclength,  Curvature,  etc.). 

2.  Given  an  a  priori  distribution,  P,  on  17,  maximize 


P(x)dx 


where  Sr  =  {a:  €  17  :  |ar  —  F|  <  c(x)}  ,  where  c(x)  is  the  radius  of  the  set 
“cut-out”  of  the  domain  by  the  path  F. 


We  note  that  this  problem  was  motivated  by  the  problem  of  computing 
optimal  search  strategies  in  the  presence  of  a  priori  knowledge  [18] .  The  function 
P  represents  any  knowledge  of  the  search  domain,  17.  Possible  choices  for  the 
optimization  would  include  minimal  arclength  and  minimal  curvature.  Also,  we 
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note  that  this  encompasses  obstacle  avoidance  when  2  is  minimized  or  the  sign 
of  P  is  changed. 


2.1  Level  Set  Formulation 


The  search  path  T  will  be  represented  by  the  0  level  set  of  a  function  <fi  : 
O  — »  R.  Given  an  initial  function  <f>(t  =  0),  and  an  energy  E(cj>)  to  be  mini¬ 
mized/maximized,  we  use  the  method  of  gradient  descent/ascent  to  arrive  at  a 
PDE  of  the  form 


_  dE 

df  ~  ~/  +  Tty’ 


(1) 


where  is  taken  from  the  Euler-Lagrange  equation.  This  PDE  is  then  evolved 
to  steady  state  resulting  in  <j>  obtaining  a  local  minimum.  Most  of  our  variational 
problems  will  be  nonconvex,  so  the  initial  choice  of  (j)  will  determine  the  local 
minimum  in  which  we  finish.  The  numerical  methods  for  solving  (1)  will  be 
discussed  later. 


2.2  Examples  of  Energies  and  PDEs 

The  energy  representing 


P(x)dx 


(2) 


will  always  be  included  in  our  variational  formulation.  First  we  assume  <j>  is  a 
weighted  signed  distance  function.  One  way  to  construct  <j>  is  to  solve 


|V0| 


1 

R(x)  ’ 


(3) 


where  R(x)  >  0,  with  boundary  condition  given  by  {(j)(x)  =  0|x  G  T}  [22], 
[19], [25], [28].  In  terms  of  the  search  problem  R  «  0  implies  that  visibility  is 
reduced  to  almost  nothing,  while  when  R  — >  oo  visibility  becomes  large.  Given 
this  assumption,  we  can  see  that  the  set 


x(Sr)  =  H(r  -  <t>)H(r  +  </>), 


(4) 


where  \  is  the  characteristic  function,  and  H  is  the  Heaviside  function.  Here  r 
is  a  constant,  and  we  assume  that  the  information  about  the  radius  of  visibility 
function  c(x)  has  been  incorporated  into  the  function  R(x). 

Thus  our  integral  (2)  can  be  written  as 


[  P(x)dx  =  f  H(r  —  <l>)H(r  +  <j>)  P{ x)  dx,  (5) 

J  Sr  ^ 

where  we  have  integrated  over  the  entire  domain  O.  To  maximize  this  integral 
we  perform  gradient  ascent  and  arrive  at  the  PDE 

4>t  =  P{x)[H (r  -  (j>)5(r  +  (j))  -  H(r  +  <t>)8{r  -  (j))},  (6) 
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where  S(y)  =  H'(y)  denotes  the  Dirac  delta  function.  Intuitively  (6)  attempts 
to  move  the  ±r  level  sets  away  from  the  0  level  set  where  P  >  0,  so  that  Jg  P 
becomes  larger  as  time  progresses. 

Another  common  energy  term  to  be  minimized  is  the  length  of  T.  We  can 
write 

|F|  =  [  5{<j>) |V0|  dx.  (7) 

■hi 

Gradient  descent  yields 

I*  =SW)V’M  =  {W)"’  (8) 

where  k  is  the  mean  curvature  of  T. 

For  certain  problems  such  as  the  tool  path  problem  one  may  want  to  control 
the  magnitude  of  k  along  the  path.  Energies  to  be  minimized  in  this  case  could 
be  of  the  form 


'K  </%(«)  dxi 


(9) 


where  g  is  a  non-negative  function  of  k  such  as  |k|p,p  >  0.  The  PDEs  resulting 
from  (9)  are  fourth  order  involving  second  partial  derivatives  of  n. 

In  general  we  will  evolve  (6)  with  its  right  hand  side  augmented  by  adding 
weighted  terms  that  are  found  from  the  other  energy  minimizations/maximizations 
that  each  particular  problem  demands.  An  example  would  be  adding  a  term 
A 8{<t>)n  to  the  right  hand  side  of  (6),  where  A  >  0  can  be  thought  of  as  a  Lagrange 
multiplier. 


3  Numerical  Methods 

3.1  Advancement  of  Time  Dependent  PDEs 

The  PDEs  found  in  section  2.2  are  generally  Hamilton- Jacobi  equations.  To  dis¬ 
cretize  them  we  construct  a  uniform  rectangular  grid  on  Q.  Viscosity  solutions 
for  these  types  of  equations  have  been  studied  well  [6],  [10]  and  numerical  meth¬ 
ods  that  converge  to  the  viscosity  solution  have  been  implemented  [7] ,  [20] .  We 
use  these  methods  to  solve  our  equations.  In  general  they  consist  of  upwind  type 
spatial  discretizations  and  explicit  Runge-Kutta  time  discretizations.  On  dfl, 
which  is  the  boundary  of  fl ,  we  use  Neumann  boundary  conditions:  dcj)/dn  =  0. 

Most  level  set  variational  problem  formulations  involve  only  one  level  set 
of  interest,  usually  the  set  To  =  {x\cl)(x)  =  0}.  Thus  the  PDEs  to  be  evolved 
only  involve  S  functions  that  have  support  localized  near  To.  However,  for 
our  problem,  (6)  includes  S  functions  whose  support  lies  near  the  sets  r±r  = 
{x\<j)(x)  =  ±r}  when  we  use  a  numerical  approximation  to  d.  When  (6)  is 
combined  with  a  Lagrange  multiplier  term  derived  from  (8),  then  we  have  a 
PDE  with  5  functions  localized  near  the  places  where  <j>  =  —  r,  0,r.  Therefore 
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three  separate  level  sets  of  <f>  must  be  evolved,  while  </>  maintains  the  criterion 
that  it  is  a  weighted  signed  distance  function  satisfying  (3)  with  To  as  the 
boundary  condition. 

If  we  try  to  simultaneously  evolve  </>  in  the  different  regions  near  ro,r±r, 
then  we  see  that  (3)  will  not  be  satisfied.  This  cannot  be  quickly  repaired  by 
reinitialization  of  cj>  to  a  signed  distance  function  because  it  is  unclear  which 
boundary  conditions  should  be  used  (i.e.  To,r_r,  or  Pr)  Interesting  motion  for 
our  variational  problem  has  occurs  near  all  three  level  sets,  and  reinitialization 
with  only  one  of  the  P2  as  the  boundary  condition  would  disregard  motion  near 
the  other  two  Tw^z. 

To  remedy  this  problem  we  propose  to  evolve  <j>  near  one  of  the  P2,  z  = 
— r,  0,  r  for  a  short  period  of  time,  then  if  z  =  ±r  we  perform  a  pseudo¬ 
reinitialization  that  approximates  the  solution  to  (3)  with  P2  as  the  boundary 
condition  (we  will  explain  in  a  moment  what  we  mean  by  ‘pseudo-reinitialization’ 
and  why  it  is  necessary),  and  if  z  =  0  we  perform  the  usual  reinitialization  to 
a  weighted  signed  distance  function  solving  (3).  Once  one  of  the  T2  has  been 
advanced  and  <j>  has  been  reinitialized,  we  move  to  another  P2  and  repeat  the 
process.  So  prior  to  each  advancement  of  a  particular  T2  we  are  forcing  <j>  to 
approximately  or  exactly  satisfy  (3). 

In  order  to  avoid  the  discretization  of  the  <5  function  we  modify  the  PDEs 
by  replacing  the  S(<f>)  with  |V</>|.  This  allows  all  level  sets  to  move  with  the 
same  speed.  However,  we  would  still  like  to  maintain  the  local  nature  of  the 
support  for  the  5  function,  so  we  also  multiply  |V^|  by  a  cutoff  function,  p(</>), 
as  described  in  [21]  with  support  over  points  where  \<j>\  <  r. 

3.1.1  Pseudo-reinitialization 

In  the  method  described  above,  if  we  were  to  use  true  weighted  signed  distance 
reinitialization  instead  of  the  pseudo-reinitialization  procedure,  we  would  en¬ 
counter  problems  where  the  reinitialization  with  boundary  condition  Ty  would 
pick  the  viscosity  solution  to  (3)  and  we  would  lose  interesting  parts  of  (f)  near 
Wi/- 

An  example  of  this  phenomenon  is  shown  in  figure  1.  Note  how  on  the  right 
plot  To  has  not  changed  location,  but  Pr  has  been  smoothed  as  the  viscosity 
solution  to  (3)  has  been  found,  eliminating  the  part  of  Pr  containing  a  cusp. 
In  both  plots  To  is  a  level  set  of  <f>  that  satisfies  (3)  with  boundary  condition 
given  by  Pr  (which  is  a  condition  that  <j>  must  satisfy),  but  one  can  see  that 
the  area  between  the  two  curves  in  much  larger  in  the  left  plot,  and  this  area 
could  be  an  energy  which  we  would  like  to  maximize.  Thus  we  would  like  to  be 
able  to  maintain  cusped  areas  of  the  T2  while  also  requiring  that  (3)  holds  with 
boundary  condition  given  by  To- 

The  method  we  use  to  achieve  this  is  denoted  pseudo-reinitialization.  For 
notational  purposes  we  will  assume  To  is  the  boundary  condition  for  reinitial¬ 
ization,  as  using  other  T2  for  the  boundary  only  require  shifts  in  the  signum 
function,  S.  The  idea  is  that  instead  of  using  the  PDE  method  of  reinitializa- 
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Figure  1:  Left:  Initial  contour  plot  of  (j>.  Dashed  line:  To.  Solid  line:  Tr.  Right: 
Contour  plot  of  (f)  after  true  reinitialization  with  To  used  as  boundary  condition, 
i.e.  solution  of  (3). 


tion  to  a  weighted  signed  distance  function  [26]  which  solves  the  equation 

*  +  =°-  m 

we  instead  solve 

<t>t  +  S(tj>)  (v<f>  ■  j^j-  —  =  0,  (11) 

where  rj  is  a  static  vector  field  found  by  taking  r)  =  \7(f>  prior  to  starting  the 
pseudo-reinitialization. 

The  difference  between  (10)  and  (11)  is  that  we  predefine  the  characteristic 
flow  directions,  t]/\t}\,  instead  of  letting  them  evolve  during  the  reinitialization, 
and  we  also  use  \rj\  instead  of  l/R(x)  to  determine  the  growth  of  </>  along  the 
characteristics.  As  long  as  <f>  has  not  changed  too  much  during  an  evolution  step 
near  To,  then  both  r]/\ri\  and  |r/|  will  be  close  to  the  quantities  they  represent  in 
(10),  which  are  \7 (j) / \\7 (f)\  and  1/R(x),  respectively. 

The  idea  behind  (11)  is  that  instead  of  allowing  the  characteristics  to  fol¬ 
low  a  rarefaction  that  would  be  found  when  taking  the  viscosity  solution  to 
(10),  as  in  figure  2,  we  instead  have  them  follow  paths  that  they  would  take  if 
reinitialization  was  being  done  with  Tr  as  the  boundary  condition  and  we  were 
finding  the  viscosity  solution  to  (10)  in  that  case.  Figure  3  shows  an  example 
of  the  paths  the  characteristics  might  take  if  Tr  had  a  protruding  bump.  These 
paths  are  quite  different  than  those  in  the  viscosity  solution  to  2,  and  could 
be  thought  of  as  being  the  characteristic  paths  of  a  nonphysical  shock  in  the 
solution  to  (10). 

To  solve  (11)  we  first  choose  tj.  This  is  done  in  a  Godunov  type  upwind 
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Figure  2:  Characteristics  of  the  viscosity  solution  to  reinitialization. 


Figure  3:  Characteristics  in  pseudo-reinitialization. 


Figure  4:  Left:  Two  contour  plots  of  <j>  at  different  times.  Solid  lines  are  initial 
contours:  To  (leftmost  curve  running  from  top  to  bottom),  Fr  (rightmost  curve 
running  from  top  to  bottom).  Dashed  lines  are  contours  of  </>  after  To  has  been 
evolved  a  short  period  of  time  (it  has  moved  slightly  to  the  left),  and  then  true 
reinitialization  has  been  performed  with  To  used  as  boundary  condition.  Right: 
Similar  contour  plots  of  <f>,  except  that  the  dashed  lines  now  have  been  calculated 
using  pseudo-reinitialization  with  Tq  as  boundary  condition. 


manner.  For  each  gridpoint  Xij  we  make  the  following  choices: 

r/i  =  maxmod(max(D~^>i!j,0),min(D'^^>itj,0))/dx, 

r]2  =  rnaxrnod(max(D~(f>itj,0),mm(Dy(f>itj,0))/dy,  (12) 

if  S(4>ij)  >  0,  and 

?7i  =  maxmod(mm(D~ (j>ij ,  0),  max(D^<^jj,  0))/dx, 

772  =  maxmod(min(D~(j)ij,Q),max.(D^(j)ij,0))/dy,  (13) 

if  <  0,  where 


maxmod{x ,  y) 


x  if  \x\  >  \y\, 
y  otherwise. 


Here  Dffaj  =  ±(<t>i± ij  -  </>;j),  and  -  <pij).  We  note  that 

more  accurate  W/ENO  methods  can  also  be  used  to  construct  r/. 

Once  7/  has  been  chosen  we  can  advance  (11)  using  upwinding  or  W/ENO 
methods  [20], [13]  in  space  and  TVD  Runge-Kutta  methods  [20]  in  time  as  it  is 
basically  a  linear  advection  equation  with  a  forcing  term. 

Figure  4  shows  an  example  of  how  pseudo-reinitialization  retains  cusped 
structures.  One  can  imagine  more  drastic  cases  where  the  point  of  the  cusped 
region  expands  to  form  a  larger  inlet  region  which  would  be  removed  during 
typical  reinitialization. 
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3.2  Topology  Preservation 

For  certain  problems  such  as  searching  it  makes  physical  sense  that  the  search 
path  To  enters  17  from  one  point  a  £  <917,  and  leaves  through  one  point  b  £  <917 
while  never  changing  topology. 

When  using  level  set  methods,  one  of  the  advantages  over  particle  tracking 
methods  is  the  ability  to  change  topology  without  user  intervention  or  adjust¬ 
ments  to  the  algorithms  being  used.  However,  in  our  case  we  may  wish  to  avoid 
these  topological  changes.  If  topology  preservation  is  required  for  a  particular 
problem  it  is  implemented  in  the  way  outlined  in  [12]. 

This  method  of  topology  preservation  in  applicable  only  to  uniform  rectan¬ 
gular  grids  and  relies  on  the  well  studied  field  of  digital  topology.  Its  implemen¬ 
tation  is  local  in  nature  and  thus  does  not  add  significant  complexity  to  the  run 
time  of  the  algorithm.  It  consists  of  a  check  of  the  signs  of  all  the  neighbors  Xk,i 
of  a  point  Xij  with  max(|fc  —  *],]/  —  j\)  <  1,  along  with  their  connectivities,  after 
all  </>itj  have  been  advanced  a  timestep.  If  a  topological  change  would  occur  in 
Fo  at  Xij  if  (f>itj  was  allowed  to  change  sign,  then  this  sign  change  is  prohibited. 
We  refer  the  reader  to  [12]  for  more  details  on  implementation  and  references 
concerning  this  technique. 

We  also  note  that  in  order  to  keep  the  points  where  F0  intersects  the  bound¬ 
ary  fixed,  we  imposed  Dirichlet  boundary  conditions  0  =  0  there.  If  these  points 
do  not  lie  on  the  uniform  grid  then  we  can  modify  the  grid  slightly  near  them 
so  that  they  are  included  in  the  discretization  of  12.  If  this  is  done  then  a  local 
method  for  advancing  the  solution  on  an  unstructured  grid  could  be  used  near 
the  points. 

3.3  Outline  of  Evolution  Procedure 

In  this  section  we  outline  the  evolution  procedure.  We  give  a  listing  of  the  steps 
taken  during  one  iteration.  The  evolution  procedure  is  repeated  until  steady 
state  is  reached. 

In  certain  steps  of  the  computational  procedure  we  shift  <j>  by  adding  or  sub¬ 
tracting  r  at  all  points  so  that  equations  such  as  (11)  involving  signum  functions 
can  be  solved  similarly  no  matter  when  they  are  called.  This  is  explained  in  this 
way  to  emphasize  that  coding  can  be  done  using  a  smaller  number  of  functions 
that  do  identical  jobs  on  shifted  versions  of  the  data.  When  this  is  done  we 
denote  the  shifted  version  of  (f)  as  </>  ±  r.  It  is  assumed  that  after  the  step  in 
question  is  completed  <fi  is  then  shifted  back  the  opposite  way  by  Sfr. 

The  evolution  loop  advancing  the  solution  from  time  t\  to  t\  +  dt  is  given 
below.  We  illustrate  the  steps  with  an  example  PDE  of  the  form: 

c At  =  P{x)[H{r  -  <f>)5(r  +  </>)-  H(r  +  <j>)8{r  -  </>)]  +  ^kS(4>),  (14) 

using  the  substitutions  of  8((j>)  functions  by  | V0|/o(0),  which  is  how  they  are 
implemented  numerically. 

1.  Find  T]  based  on  <j>  +  r  using  (12),  (13). 
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2.  Evolve  from  time  t.\  to  t\  +  dt  the  points  near  where  </>  =  —  r,  i.e.  evolve  all 
PDE  terms  with  <5(r  +  </>)  in  them,  e.g.  4>t  =  P(x)  H{r  —  4>)\V(f>\p(r  +  (f>). 

3.  Pseudo-reinitialize  using  (11)  with  (f>  +  r. 

4.  Find  77  based  on  </>  —  ?’  using  (12),  (13). 

5.  Evolve  from  time  t\  to  t\  +  dt  the  points  near  where  </>  =  r,  i.e.  evolve  all 
PDE  terms  with  8{r  —  (f> )  in  them,  e.g.  (j>t.  =  P(x)  H{r  +  <fi)\V(f>\p(i'  —  <fi). 

6.  Pseudo-reinitialize  using  (11)  with  <f>  —  r. 

7.  Evolve  from  time  t\  to  t\  +  dt  the  points  near  where  cf>  =  0,  i.e.  evolve  all 
PDE  terms  with  the  <5(</)  in  them,  e.g.  </t  =  p,K\Vcj)\p  ((/)). 

8.  Reinitialize  to  a  weighted  signed  distance  function  from  the  0  level  set 
using  (10). 

9.  Go  to  step  1. 


4  Numerical  Simulations 


In  this  section  we  present  some  numerical  simulations.  The  PDE  we  evolve  to 
steady  state  is 


4>t  =  P(x)[H(r  -  (j))5{r  +  </>)  -  H(r  +  cj))S(r  -  </>)]  +  pn5((j>),  (15) 

using  the  methods  mentioned  above.  The  domain  £T  is  [ —  1 , 1] 2  for  all  problems, 
discretized  in  a  uniform  rectangular  grid  with  dx  =  dy  —  1/80.  As  mentioned 
above  Neumann  BCs  d(f>/dn  =  0,  are  used.  Topology  preservation  is  used  in 
all  figures  except  for  figure  12.  After  replacing  the  S  functions  with  \\7(j)\p((j))  a 
conservative  estimate  on  the  CFL  condition  for  the  problem  is 


dt  max 


M 

dx 


\p\  _ 

dy  1  dx2 


(16) 


where  we  use  the  max  applied  to  the  P  and  n  terms  individually  instead  of  to 
their  sum  because  we  are  splitting  the  evolution  procedure. 

For  the  individual  examples  we  do  not  explicitly  write  the  initial  conditions, 
but  rather  shown  them  in  contour  plots.  The  way  they  are  constructed  is  by 
determining  an  initial  contour  To(t  =  0)  and  finding  an  arbitrary  function  that 
takes  r0(f  =  0)  as  its  0  level  set,  and  then  running  reinitialization  over  the  entire 
domain  fi  with  To(t  =  0)  as  the  boundary  condition. 

For  each  example  we  show  the  initial  and  steady  state  solutions  as  well  as 
a  plot  of  the  energy  over  time.  This  energy  function,  defined  by  taking  (5)  —  p 
times  (7),  is  what  we  hope  to  maximize.  The  6  and  Heaviside  functions  used  are 
the  compactly  supported  ones  given  in  [4],  with  support  parameter  e  =  2 dx.  It 
should  be  noted  that  a  more  accurate  numerical  construction  of  these  singular 
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Figure  5:  Initial,  steady  state  and  energy  vs.  timesteps  plots.  Constant  visibility 
function  R(x). 


Figure  6:  Initial,  steady  state  and  energy  vs.  timesteps  plots.  Constant  visibility 
function  R{x). 


functions  can  be  found  in  [9]  should  a  more  exact  measure  of  the  energy  be 
needed. 

In  all  examples  the  path  boundary  level  sets  are  r  =  ±0.25. 

In  figure  5  we  show  the  initial,  steady  state  and  energy  plots  for  an  example 
with  constant  R(x)  =  1.  The  circular  objects  have  P(x)  =  40,  while  P(x)  =  0 
otherwise.  The  regularization  coefficient  y,  =  0.5. 

In  figure  6  we  show  the  initial,  steady  state  and  energy  plots  for  an  example 
with  constant  R(x)  =  1,  but  this  time  we  use  a  different  initial  condition.  The 
circular  objects  have  P(x)  =  40,  while  P(x )  =  0  otherwise.  The  regularization 
coefficient  //  =  1 . 

In  figure  7  we  show  the  initial,  steady  state  and  energy  plots  for  an  example 
with  varying  R(x).  In  the  lower  half  of  the  plane  R(x)  =  0.6,  and  in  the  upper 
half,  R(x)  =  1.  The  circular  objects  have  P(x)  =  40,  while  P(x)  =  0  otherwise. 
The  regularization  coefficient  /i  =  0.5. 

In  figure  8  we  show  the  initial,  steady  state  and  energy  plots  for  an  example 
where  P(x)  is  nonzero  in  a  larger  region  of  interest.  Within  a  circle  of  radius 
0.9  we  set  P{x)  =  20  +  20  G(x,y),  where  G{x,y )  is  a  compactly  supported 
approximation  to  a  S  function  with  radius  of  support  0.9,  see  figure  9.  In  figure 
8  the  boundary  of  the  region  with  nonzero  support  is  shown.  In  this  example 
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Figure  7:  Initial,  steady  state  and  energy  vs.  timesteps  plots.  Varying  visibility 
function  R(x). 


Figure  8:  Initial,  steady  state  and  energy  vs.  timesteps  plots.  Gaussian  type 
P{x). 


R{x)  =  0.5  in  the  entire  domain,  and  /i  =  0.5.  In  this  example  the  steady  state 
solution  is  the  simple  back  and  forth  sweeping  pattern. 

In  figure  10  we  show  the  initial,  steady  state  and  energy  plots  for  an  example 
where  object  avoidance  is  the  objective.  Here  P(x)  =  —40  in  the  circular  objects 
and  0  elsewhere.  In  this  example  R(x)  =  0.5  in  the  entire  domain,  and  /j  =  1. 

In  figure  11  we  show  what  can  happen  if  the  initial  conditions  are  too  close  to 
a  local  minimum  that  is  far  form  the  global  minimum.  The  problem  parameters 
are  identical  to  those  of  figure  10,  but  a  different  initial  condition  is  imposed. 

In  figure  12  we  show  the  initial,  intermediate,  steady  state  and  energy  plots 
for  an  example  where  topology  change  is  allowed.  Here  P(x)  =  40  in  the  circular 
objects  and  P(x )  =  —40  elsewhere.  In  this  example  R{x)  =  0.5  in  the  entire 
domain,  and  /r  =  0.02. 

Note  how  concentric  ellipses  have  formed  within  each  of  the  elliptical  regions 
where  P  >  0. 
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Figure  9:  P{x)  used  for  figure  8. 


Figure  10:  Initial,  steady  state  and  energy  vs.  timesteps  plots  for  object  avoid 
ance.  Solution  converges  to  global  minimum. 


Figure  11:  Initial,  steady  state  and  energy  vs.  timesteps  plots  for  object  avoid 
ance.  Solution  converges  to  local  minimum. 
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Figure  12:  Top  row:  initial,  and  intermediate  state  plots.  Bottom  row:  steady 
state,  and  energy  vs.  timesteps  plots.  Topology  change  is  allowed  in  this  exam¬ 
ple. 
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5  Conclusion 


We  have  presented  a  level  set  based  algorithm  for  solving  a  variational  approach 
to  path  planning.  Some  key  features  of  this  algorithm  are  the  energy  integrals 
used  to  define  the  search  criteria,  the  splitting  technique  used  to  advance  the 
PDEs,  and  the  pseudo-reinitialization.  The  number  of  timesteps  needed  to 
converge  to  steady  state  are  on  the  order  of  hundreds,  and  for  2  dimensional 
problems  the  runtime  can  be  made  close  to  real  time  for  applications. 

The  energy  integrals  used  are  very  basic  and  encompass  general  properties 
that  are  desirable  in  many  path  planning  problems.  However,  they  are  not  ex¬ 
haustive  and  more  complicated  energies  based  on  functionals  of  curvature  and 
other  path  properties  can  be  constructed.  Solving  some  of  these  types  of  varia¬ 
tional  formulations  would  be  straightforward,  but  there  are  definite  challenges 
ahead  in  this  area. 

An  analytic  proof  of  the  numerical  convergence  of  the  splitting  and  pseudo¬ 
reinitialization  techniques  to  a  true  maximization  of  the  energy  is  also  needed. 
While  these  techniques  give  quantitative  results  that  indicate  energy  maximiza¬ 
tion,  it  would  be  useful  to  have  a  better  understanding  of  their  numerical  prop¬ 
erties. 

Extensions  to  3  dimensions  would  also  be  useful,  but  in  3d  the  path  To  is  a 
codimension-2  set  and  therefore  either  requires  the  use  of  2  level  set  functions, 
or  a  new  way  of  being  tracked.  Also,  the  splitting  scheme  in  3d  may  become 
more  complicated. 

Some  other  problems  which  we  have  not  approached  but  are  feasible  for 
future  research  are:  multiple  non-intersecting  paths,  time  dependent  parameters 
such  as  R,  P,  /i,  paths  passing  through  multiple  prescribed  points,  and  self- 
intersecting  paths.  Level  set  motion  with  interfaces  having  positive  width  may 
also  have  applications  in  image  processing,  such  as  segmentation  of  thin  objects, 
and  also  in  areas  of  physics  where  the  dynamics  may  demand  that  an  interfacial 
width  take  nonzero  measure. 
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